/*------------------------------------------------------------------------
 *  Copyright (C) 2011  Luis M. de la Cruz Salas
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

#ifndef _TWOPHASEEQUATION_H_
#define _TWOPHASEEQUATION_H_

#include "Tuna.hpp"
#include "Equations/GeneralEquation.hpp"

namespace Tuna {

template <typename T>  class TwoPhaseEquation;
  
template <typename T>  
struct Typeinfo <TwoPhaseEquation<T> > {
  typedef typename Typeinfo<T>::prec_t prec_t;
  enum { Dim = Typeinfo<T>::Dim };
};

/*! 
 ***************************************************************************
 *  This class derives from GeneralEquation through B & N trick.
 *  In this class some particular variables and methods for scalar PDE's
 *  are defined.
 ***************************************************************************
 *  \author  Luis M. de la Cruz [ Mon May  4 16:20:22 CDT 2009 ]
 ***************************************************************************
 */

template<typename Tscheme>
class TwoPhaseEquation : public GeneralEquation<TwoPhaseEquation<Tscheme> >
{
  using GeneralEquation<TwoPhaseEquation<Tscheme> >::name;

public:
  typedef typename Typeinfo<Tscheme>::prec_t prec_t;
  typedef typename TunaArray<prec_t, Typeinfo<Tscheme>::Dim >::tiny floatTinyArray_t; 
  typedef typename TunaArray<int, Typeinfo<Tscheme>::Dim >::tiny intTinyArray_t;
  typedef typename TunaArray<prec_t, Typeinfo<Tscheme>::Dim >::huge ScalarField;
  typedef SparseMatrix< Diagonal< prec_t, Typeinfo<Tscheme>::Dim > > DiagMat;


  TwoPhaseEquation() : GeneralEquation<TwoPhaseEquation<Tscheme> >() { }

/*!
 *  Main constructor.
 *  The scalarfield \c phi_0 defined here is a reference to a global field 
 *  variable. Here \c phi_0 is used as an old value of that field variable. 
 */
  TwoPhaseEquation(ScalarField &phi_global, 
		 DiagMat &matrix, ScalarField &b,
		 const floatTinyArray_t& deltas)
    : GeneralEquation< TwoPhaseEquation<Tscheme> >(phi_global, matrix, b, deltas)
  {
    Dim = Typeinfo<Tscheme>::Dim;
    phi_0.reference(phi_global);
    name = "TwoPhaseEquation";
    mu_o = -1;
    theta = 1;
  }
  ~TwoPhaseEquation() { }; 
  
/*!
 *  The Curiously Recursive Template Pattern (CRTP) is used.
 *  asDerived() gets an returns a reference to the derived class.
 */
  inline Tscheme& asDerived() { return static_cast<Tscheme&>(*this); }
  
//  
// Set functions used to initialize the parameters.
//
  inline void setPermeability(prec_t permeability) { k = permeability; }
  inline void setSrw(prec_t srw) { Srw = srw; }
  inline void setSro(prec_t sro) { Sro = sro; }
  inline void setInjection(prec_t i) { injection = i; }
  inline void setViscosity_w(prec_t vw) { mu_w = vw; }
  inline void setViscosity_o(prec_t vo) { mu_o = vo; }
  inline void setPorosity(prec_t po) { porosity = po;}
  inline void setTheta(prec_t th) { theta = th;}
  inline void setPthreshold(prec_t p_t) {p_threshold = p_t;}

  inline void setSaturation(const ScalarField& ss) { S.reference(ss); }
  inline void setPressure(const ScalarField& pp) { p.reference(pp); }

//  
// Get functions
//
  inline ScalarField& get_Phi_0() { return static_cast<ScalarField&>(phi_0);}
  inline prec_t getPorosity() { return porosity; }

/*!
 *  Call to the functions defined in the derived classes.
 */
  inline bool calcCoefficients()
  {
    if(Dim == 1)
      return asDerived().calcCoefficients1D();
    else if (Dim == 2)
      return asDerived().calcCoefficients2D();
    else if (Dim == 3)
      return asDerived().calcCoefficients3D();
  }

  inline bool updateBoundaries()
  {
    if(Dim == 1)
      return asDerived().updateBoundaries1D();
    else if (Dim == 2)
      return asDerived().updateBoundaries2D();
    else if (Dim == 3)
      return asDerived().updateBoundaries3D();
  }

  inline void printInfo() {
    std::cout << "\n +-----+"
	      << "\n | TwoPhaseEquation< ";
    asDerived().printInfo();
    std::cout << "< " << Dim << "D  > >"
	      << "\n +-----+"
	      << "\n | Permeability (k)= " << k 
	      << "\n | Porosity = " << porosity	     
	      << "\n | Water Residual Sat (Srw) = " << Srw
	      << "\n | Oil Residual Sat (Sro) = " << Sro
	      << "\n | Water viscosity (mu_w) = " << mu_w
	      << "\n | Oil viscosity (mu_o) = " << mu_o
	      << "\n | Pore size distrib (theta) = " << theta
	      << "\n | Threshold pressure (Pt) = " << p_threshold
	      << "\n +-----+";  
  }

  inline void setRank(int r) { rank = r; }
  inline void setSize(int s) { size = s; }
  
protected:
  int Dim, rank, size;
  prec_t k, Srw, Sro, mu_w, mu_o, porosity, theta, p_threshold;
  prec_t injection, extraction;

public:

  ScalarField S, p;
  ScalarField phi_0;
};


}  // namespace Tuna

#endif //_TWOPHASEEQUATION_




